#Compare means of the four items by treatment
#Plot the Means
library(Rmisc)
library(DescTools)
library(ggplot2)
study2$treatment2<-as.factor(study2$treatment)
table(study2$treatment2)
#drop the cases with missing values of the authoritariansim measures.

study2<-subset(study2, Obedience>0 & Considerate>0 & Curiosity>0 & Independence>0)

plot(study2$treatment2)

levels_temp=c("Boy", "Control", "Girl")


##Plot the means by condition
ind_mean <- summarySE(study2, measurevar="Independence",
                      groupvars="treatment2",
                      na.rm=TRUE, conf.interval=.95)
ind_mean<-ind_mean[-4,]

ind_mean$test<-c("Boy", "Control", "Girl")
##recode to zero one
ind_mean$Ind2<-ind_mean$Independence-1
summary(study2$Independence)

elder_plot<-ggplot(ind_mean, aes(x=factor(test, level =c('Control','Boy', 'Girl')), y=Ind2, group=1)) +
  geom_line() +
  geom_errorbar(width=.1, aes(ymin=Ind2-ci, ymax=Ind2+ci), colour="red", data=ind_mean) +
  geom_point(shape=21, size=1, fill="white")+
  xlab("Treatment")+
  ylab("Proportion choosing \n'Respect for Elders'")+
  ylim(0,1)+theme_bw()

elder_plot
rep.Independence.aov<-aov(Independence~as.factor(treatment2), data=study2)
summary(rep.Independence.aov)
TukeyHSD(rep.Independence.aov)



obed_mean <- summarySE(study2, measurevar="Obedience", groupvars="treatment2",
                       na.rm=TRUE, conf.interval=.95)

obed_mean<-obed_mean[-4,]
##recode to zero one
obed_mean$obed2<-2-obed_mean$Obedience
obed_mean$test<-c("Boy", "Control", "Girl")
obed_mean
obed_plot<-ggplot(obed_mean, aes(x=factor(test, level =c('Control','Boy', 'Girl' )), y=obed2, group=1)) +
  geom_line() +
  geom_errorbar(width=.1, aes(ymin=obed2-ci, ymax=obed2+ci), colour="red",data=obed_mean) +
  geom_point(shape=21, size=1, fill="white") +
  xlab("Treatment")+
  ylab("Proportion choosing \n'Obedience'")+
  ylim(0,1)+theme_bw()
rep.Obedience.aov<-aov(Obedience~as.factor(treatment), data=study2)
summary(rep.Obedience.aov)
TukeyHSD(rep.Obedience.aov)


cur_mean <- summarySE(study2, measurevar="Curiosity", groupvars="treatment2",
                      na.rm=TRUE, conf.interval=.95)

cur_mean
cur_mean<-cur_mean[-4,]
##recode to zero one
cur_mean$cur2<-cur_mean$Curiosity-1
cur_mean$test<-c("Boy", "Control", "Girl")

cur_plot<-ggplot(cur_mean, aes(x=factor(test, level =c('Control','Boy', 'Girl')), y=cur2, group=1)) +
  geom_line() +
  geom_errorbar(width=.1, aes(ymin=cur2-ci, ymax=cur2+ci), colour="red", data=cur_mean) +
  geom_point(shape=21, size=1, fill="white") +
  xlab("Treatment")+
  ylab("Proportion choosing \n'Good Manners'")+ylim(0,1)+theme_bw()

rep.Curiosity.aov<-aov(Curiosity~as.factor(treatment), data=study2)
summary(rep.Curiosity.aov)
TukeyHSD(rep.Curiosity.aov)


con_mean <- summarySE(study2, measurevar="Considerate", groupvars="treatment2",
                      na.rm=TRUE, conf.interval=.95)

con_mean
con_mean<-con_mean[-4,]
##recode to zero one
con_mean$con2<-con_mean$Considerate-1
con_mean$test<-c("Boy", "Control", "Girl")

wb_plot<-ggplot(con_mean, aes(x=factor(test, level =c('Control','Boy', 'Girl')), y=con2, group=1)) +
  geom_line() +
  geom_errorbar(width=.1, aes(ymin=con2-ci, ymax=con2+ci), colour="red", data=con_mean) +
  geom_point(shape=21, size=1, fill="white") +
  xlab("Treatment")+
  ylab("Proportion choosing \n'Well-Behaved'")+ylim(0,1)+theme_bw()

rep.Considerate.aov<-aov(Considerate~as.factor(treatment), data=study2)
summary(rep.Considerate.aov)
TukeyHSD(rep.Considerate.aov)

library(cowplot)
pdf(file="Figure2.pdf")
plot_grid(elder_plot, obed_plot, cur_plot, wb_plot)
dev.off()

Author_mean <- summarySE(study2, measurevar="auth", groupvars="treatment2",
                         na.rm=TRUE, conf.interval=.95)

Author_mean$Author<-Author_mean$auth
Author_mean$test<-c("Control", "Boy", "Girl")

plot_Author<-ggplot(Author_mean, aes(x=factor(test, level =c('Control','Boy', 'Girl')), 
                                     y=Author, group=1)) +
  geom_line() +
  geom_errorbar(width=.1, aes(ymin=Author-ci, 
                              ymax=Author+ci), colour="red", 
                data=Author_mean) +
  geom_point(shape=21, size=1, fill="white") +
  xlab("Treatment")+
  ylab("Authoritarianism measure")+ylim(0,4)+ theme_bw()
pdf(file="Figure5.pdf")
plot_Author
dev.off()
Author.aov<-aov(auth~as.factor(treatment), data=study2)
summary(Author.aov)
TukeyHSD(Author.aov)



